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ABSTRACT 



An algorithm is developed which enables the computation 
of horizontal range and/or depth from a submerged sound 
source, using ray acoustics and the Lloyd's mirror 
interference effect. The solution is based on Snell's law 
and involves integrating multipath sound rays to find the 
difference in length between the direct and surface 
reflected sound paths from the source to the receiver. This 
difference in path length is directly related to the 
observed Lloyd's mirror interference pattern. 

No assumptions as to the mathematical characteristics of 
the sound velocity profile (SVP) are made nor are far field 
approximations used. The solution is as accurate as the SVP 
data input to the problem. 

A general computer flow chart and basic language program 
are provided to allow local commands the capability of pre- 
mission planning based on specific operating area 
environmental information. 
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I. 



INTRODUCTION 



Sound propagation between a source and a receiver in the 
ocean may occur by many different ray paths, each dependent 
on existing physical and acoustical conditions. Crucial to 
the determination of which sound paths are of importance in 
any given problem are factors such as the sound velocity 
profile, the depth of the water, the bottom contour and 
composition, the sea surface conditions and the intensity 
and frequency of the sound. 

This thesis deals with two specific paths of sound 
propagation from a source to a receiver. The first, called 
the direct path, does not touch the surface or bottom of the 
sea as it travels from the source to the receiver. The 
second path reflects off the sea surface before reaching the 
receiver and is called the reflected path. 

The following explanation of the reception of sound from 
a single source via these two distinct paths is taken 
directly from Urick; [Ref. 1: pp. 120-123] 

"If the sea surface were perfectly smooth, it would 
form an almost perfect reflector of sound. The intensity 
of sound reflected from the smooth sea surface would be 
very nearly equal to that incident upon it." 

He continues: 

"A criterion for the roughness or smoothness of a 
surface is given by the Rayleigh parameter, defined as 
R= kH sin0 ;where k is the wave number, (27 t)/a ( A is the 
wave length of the sound incident on the surface), H is 
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the rms (sea) wave height (crest to trough), and 0 is the 
grazing angle". 




Figure 1-1. Rayleigh Parameter (R) for Describing Sea 

Surface Roughness. 

"When R<<1, the surface is primarily a reflector and 
produces a coherent reflection at the specular angle equal 
to the angle of incidence. When R>>1, the surface acts 
as a scatterer, sending incoherent energy in all 
directions. When certain theoretical assumptions are 
made, (See: P. Beckmann, and A. Spizzichino: Scattering of 
Electromagnetic Waves From Rough Surfaces , p. 93, The 
Macmillan Company, N.Y., 1963.) the (pressure) amplitude 
reflection coefficient m of an irregular surface (defined 
as the ratio of the reflected or coherent pressure 
amplitude of the return to the incoherent pressure 
amplitude) can be shown to be given by the simple 
expression m = exp (-R)." 

Urick goes on to explain the image interference. 

"When the sea surface is not too rough, it creates an 
interference pattern in the underwater sound field. This 
pattern is caused by constructive and destructive 
interference between the direct and surface reflected 
sound and is called the Lloyd mirror, or image 
interference, effect." 

In Urick's sound propagation in the sea [Ref. 2: p. 9-5], he 

continues the discussion of surface reflection interference: 

"This sound field may be divided into three parts, 1) a 
NEAR FIELD close to the source in which the image source 
is too far away, and the reflected sound is too weak to 
produce appreciable interference, 2) an INTERFERENCE FIELD 
in which there are strong loops and nulls in the signal 
received by a receiver moving outward in range, and 3) a 
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FAR FIELD in which there is an increasingly out-of-phase 
condition between source and image and the intensity falls 
off as the inverse fourth power of the range (transmission 
loss increases as 40 log r)." 

The interference effects felt at the receiver is 
expressable as the ratio of the combined direct and 
reflected sound intensities to the direct sound intensity 
alone. 

Stated mathematically: 



Intensity Ratio = 

or more precisely 
Intensity Ratio = 



^ t ot al 



^direct 



^(direct'*’ reflected) 



^direct 




2m 




COS(wT ) 



Where Rq is the direct ray path length, Rr is the reflected 
path length, m is the surface pressure reflection 
coefficient, w (where w = 27rf) is the angular frequency of 
the sound, and r is the time delay caused by the difference 
in the two path lengths. A formal derivation of the 
intensity ratio is given in Appendix A. A more 
comprehensive treatment of the effects of surface roughness 
on acoustic waves is given by Clay and Medwin [Ref. 3]. 
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II. APPLICATION of LLOYD MIRROR INTERFERENCE PHENOMENON to 
LOCATE the HORIZONTAL RANGE and DEPTH of the SOUND 
SOURCE FROM the RECEIVER 

Tactical tracking information from the Lloyd's mirror 
phenomenon has been used operationally by the Navy for many 
years. In an excellent report on the analysis of Lloyd's 
mirror phenomenon under constant velocity gradient 
conditions, Shudde [Ref. 4: p. 3], discusses the difficulty 
of solving the surface reflection interference problem. 
Accurate error estimation for any acoustical condition other 
than the isovelocity case is not practical unless a large 
computer is used. Since Shudde's 1977 report, the 
advancements in memory size, and speed of operation of small 
portable desk top computers now allow an accurate solution 
of source position for all sound velocity conditions. The 
following mathematical solution is based on the lengths of 
the sound paths determined from Snell's Law. As explained 
by Urick, [Ref. 1 : p. 116], 

"One of the most important practical results of ray 
theory is Snell's Law, which describes the refraction of 
sound rays in a medium of variable velocity." 

The solution is not dependent on any isovelocity or constant 

gradient conditions other than the reasonable assumption 

that no horizontal sound speed gradients exist throughout 

the operating area. The mathematical solution presented in 

Section III makes no assumption as to the mathematical 
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characteristics of the sound velocity profile nor does it 
use far field approximations. Past solutions depended on 
specific sound velocity profiles, such as an isovelocity or 
"straight line" profile which resulted in significant 
positional errors for some tactical applications when the 
actual sound velocity profile did not match the predicted 
model . 

In applying the solution to the Lloyd's mirror 
interference effect one must understand the relationship of 
the propagating sound to the interference pattern. Coppens, 
Sanders and Dahl [Ref. 5 : p. 55 ] explained, 

"...the noise radiated by underwater vehicles consists 
of two contributions: broadband noi se (similar to the 

hissing found between FM radio stations) and tonals 
(single frequency tones, like whistles). Thus, the 
radiated noise consists of all frequencies between certain 
broad limits with certain specific frequencies occurring 
with high intensity. Each of these frequencies, whether 
of broadband noise or tonal, has its own interference 
pattern. 

Because each frequency has its own period and 
wavelength, each will have a different phase delay for 
propagation over the direct and surface reflected paths. 
As a result, the ranges between source and receiver for 
which the surface reflected signal tends to cancel the 
direct signal will be different for each different 
frequency. Thus, the ranges at which different 
frequencies will have nulls at the receiver will be 
different. This means that as the range between source 
and receiver changes, the frequencies which have nulls at 
the receiver will also change." 

Figure II-1 shows a schematic time vs. frequency plot of the 

interference pattern as the distance between the source and 

receiver is decreased and then increased. The distinct 

tonals will most readily show the effects of Doppler shift 
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(apparent frequency change) as minimum separation range, is 
passed (see t = 0 in Figure II-1). The light (unshaded) 
hyperbolic shaped regions centered about t = 0 are known as 
the Lloyd's mirror interference pattern for the broad band 
noise. The frequency values F-| , F 2 , ^2' ^4 i^®prssent the 
minimvun frequency of each hyperbola. The actual frequency 
from the sound source that is being nullified by the surface 
reflection interface varies with time as the hyperbolas are 
traced out. The measured separation in hertz between any 
two successive nulls at any time t, is termed the difference 
frequency, A f , and is the frequency value required in the 
problem solution discussed in the following Section. (See 
Section III.D.2) 



1 4 




Figure II -1 . 



Frequency vs. 
Interference 



Time Display of 
Pattern 



Lloyd 's 



Mirror 
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III. RANGE vs. DEPTH SOLUTION USING the LLOYD ' S MIRROR 

INTERFERENCE PHENOMENON 



A. THEORETICAL BACKGROUND 
1 . Snell ' s Law 

Snell's Law, [Ref. 1; p. 116, Ref. 6: p. 120], is a 
result of Fermat's minimum time principle which states that 
each ray will travel the path requiring minimum time of 
travel from source to receiver. 

Mathematically, Snell's Law is expressed as: 



sin 0 

Snell Constant = 

c(z) 



Where 0 is the angle that the ray makes with the vertical 
and c(z) is the sound velocity for a given depth z. (See 
Figure III-1 ) 



^initial 




Figure III-1. Snell's Law Applied to Initial Sound Source 

Ray Angle. 
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The sound ray maintains the same Snell Constant value 
throughout its transmission path. Sound velocity changes 
result in an increasing or decreasing angle 0. 

2 . Ray Path Integration 

Expanding the diagram of the sound ray, the 
differential raypath length, ds, can be described as a 
function of the differential change in depth, dz, and the 
differential change in horizontal range, dr. (Figure III-2) 




Figure III-2. Pythagorean Theorem Applied to Ray Path 

Length. 

To calculate the actual path length of the sound 
ray, smaller and smaller segments of path length, ^s, 
horizontal range, Ar, and depth, Az are used until the 
differential segments ds, dr, and dz are approached. This 
tends to smooth the curve of the ray into a sequence of 
straight line segments. The path length traveled by a 
particular ray can then be determined by integrating between 
the source depth and receiver depth. The possibility of the 
direct ray crossing the same depth more than once is not 
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considered a problem in this solution due to the short 
ranges involved in Lloyd's mirror observations and the long 
ranges required for a sound ray to curve back through the 
same depth before reflecting off the surface. To handle the 
fact that the surface reflected ray does cross the same 
depth twice in the ranges pertinent to this discussion and 
to avoid errors in integration, the concept of a geometric 
equivalent image [Ref. 1: p. 123], of the source is used. 
The image reflected path and the actual reflected path are 
the same length. (Figure III-3) 




+ Range 
Meters 



Figure III-3. Actual and Image Sound Paths From Source to 

Receiver for a Given Sound Velocity Profile 
( SVP ) . 



The curved sound path shown in Figures III -3 and 
III-4 is due to the change (in general) in the sound 
velocity as the depth changes. It is possible to show that 
the path from A to B is the same length as the path from A 
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I 

to B . Also, by inspection of the sketch, the horizontal 

I 

range from A to B or A to B is the same. Therefore, for 
the reflected path calculations of path length and 
horizontal range, the image reflection depth of will be 

used . 

B. HORIZONTAL RANGE SOLUTION 
1 . Direct Path 

The horizontal range traversed by the direct path 
sound ray can be calculated from Snell's Law and the 
Pythagorean Theorem for right triangles applied to Figure 
III-2. 



From 




the Snell Constant, K-^ is: 



sin 

C(Zj^) 




This is the Snell constant for the particular direct path 
ray that ultimately travels from the source to the receiver, 
but initially starts at depth z^ and moves in direction 
Using the Pythagorean Theorem on Figure III-2 

dr 

sin 0 = 




(dz)2 



Solving 




for sin 0, 



sin 0 = K-j c(z) 



where 0 = variable for the 
particular direct path ray 
tagged with Snell constant 
Ki, 




19 



Du( 

Uoi 



setting and equal to each other yields: 



Ki c(z) 



dr 




(dz)2 




The right side of can be simplified by dividing 

numerator and denominator by dr. 



K-| c(z) = 



1 




2 



Solving for dz/dr follows in steps: 




1 + 




dz 

dr 



Kl 


c(z) 






1 




[K^ 


c(z) 


]2 


1 - [K^ 


c(z ) 


]2 


[Ki 


c(z) 




J' - ' 


Ki c( 


z)]2 



[K-| c(z)] 
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dz 



© 



dr 



yi - [Ki c(z)]2 
Ki C(z) 



Now, solving for dr, the incremental horizontal 
range yields; 



Kt c(z) 

dr = + dz (?) 

y 1 - [Kt c(z) ]2 



Integrating (^ with respect to depth, with the 
limits of integration being the depth of the source and the 



depth of the receiver, yields; 



R (direct path) 



Source 

Depth 




Actual 

Receiver 

Depth 

(along path of ray 
"tagged" ) 




The "Range" is always positive and thus the 
appropriate + or - sign is used to make this integral always 
positive regardless of the integration limits. For ease of 
discussion in the rest of this solution, the receiver will 
be assumed to be at a shallower depth. 
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2. Reflected Path 



Recalling the image depth discussion and noting a 
different angle, 6^, for the reflected ray is required 
(Figure III-4), which results in a different Snell constant, 
K 2 , where 

sin (0j^) 



It follows directly that can be rewritten for the 

reflected horizontal range as: 




R (reflected path) 



Source 



Depth 




Receiver 

Depth 



dz 
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Figure III-4. Common Direct and Reflected Sound Paths From 

Source to Receiver. 

C. SOUND RAY PATH LENGTH SOLUTION 
1 . Direct Path 

The direct path ray length can be calculated by 
returning to Figure III-2, applying the Pythagorean Theorem 
and noting: 

ds = \/(dr)^ + (dz)^ 



realizing the path length S from receiver to source is given 



by: 



Source 

Depth 

Receiver 

Depth 



Source 

Depth 

j sj (dr)2 + (dz)2 

Receiver 

Depth 
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Which can be simplified by multiplying by 



dz 



\J dz2 

as follows: 



Source 

Depth 




Receiver 

Depth 



Substituting the equivalent of 




/dr\ 

I — I from the reciprocal of 
\dz / 



S 



Source 

Depth 



(Direct Path) 




Receiver 

Depth 



+ 



^ c ( z ) ^ 

1 - c(z)2 



dz 



Source 

Depth 



1 - ^ c(z)^ + K-| ^ c(z)^ 

1 - [K-| c(z)]^ 

Receiver 

Depth 



S (Direct Path) = 



h 



dz 
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Source 

Depth 



S (Direct Path) 




1 



- [Ki c(z)]2 



dz 




Actual 

Receiver 

Depth 



2 . Reflected Path 

Using the same arguments about the reflected path as 
presented for horizontal range and substituting K 2 for K-| 
and image receiver depth into ( 11 ) yields: 



S (Reflected Path) 



Source 

Depth 




Image 

Receiver 

Depth 



dz 




D. EXPLANATION OF SOLUTION 
1 . Some Observations 

Equations and are path integrals for the 
direct ray, using the same Snell constant and integration 
limits. Further, equations and are path integrals 
of the reflected ray and thus use a different Snell constant 
from the direct path integrals, but use the same integration 
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limits respectively for the horizontal range and ray path 
lengths. 

2 . Solution 

The key to the eventual numerical solution is the 
equality between the direct path horizontal range and the 
reflected path horizontal range. Although the actual path 
length of the reflected ray is longer than the direct path 
ray between common source and receiver, the horizontal 
distance traveled by both is the same. 

By calculating all of the initial ray angles that 
are possible for each depth of the source of interest, two 
horizontal range arrays, one for the direct path and one for 
the reflected path, are formed. Each of these arrays is 
linked respectively to a direct and reflected ray path 
length array by the corresponding Snell Constant used to 
determine both the horizontal range and sound ray length for 
each path, (See Figure III -4). For example, suppose the 
initial direct path ray angle of 15 degrees from the 
vertical covered a horizontal range of 20 meters in 
traveling from the source to the receiver. Similarly, the 
initial reflected path ray angle of 8 degrees from the 
vertical covered the same horizontal range of 20 meters from 
the source to the receiver. Then the difference between the 
reflected ray path length for 8 degrees and the direct ray 
path length for 15 degrees could be computed (Figure III-5). 
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Figure III-5. Example of Reflected Path Ray Interfering 

With Direct Path Ray at the Receiver. 

It is this difference in ray path length, 
= S - S , for a wide band of frequencies that causes the 
time delayed phase shift sensed at the receiver and known as 
the Lloyd's mirror interference phenomenon. Since both the 
direct and reflected rays are traveling at the same sound 
velocity at the receiver, the ray path length difference can 
be directly converted to the frequency difference, Af as 
observed at the rece iver and seen in the Lloyd's mirror 
interference pattern. 

The intensity ratio at the receiver is given by 

/Rd\ 

2 Ml ) cos ( Ci> T ) 

WJ 



Intensity ratio = 1 + 



C-;) 
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The ratio reaches a minimum value corresponding to 
maximum interference when cos(w t) =1. This occurs when 
« T = 2 T or any integer multiple of 2 jt . 

Note that X = c/f where X = wavelength (m) 

c = phase speed (m/s) 
f = frequency in (hz) 

k = wave number {m“^ ) 

and k = 2 n- / X 

Now let A0 equal the phase difference seen at the 
receiver due to the two path length differences and defined 
as , 



A0 = k (Sj^- Sq), where 
= Path length of the reflected ray and 
Sq = Path length of the direct ray. 

Let AS = (Sjj- Sq) 

So that the first null occurs when 



A0 = k(AS)=2x=wr 

Sok (AS) = WT = 2 r = 2ir{AS)/X or 
A S = 2 tX/2 t = X 
therefore AS = c/f-j 

so the first null occurs at frequency 
fl = c/ AS 



The second null occurs at 
w T = 4 r so that 
AS = 4 ir \ 1 2 IT = 2 X = 2c/f2 

therefore second null occurs at frequency 
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±2 = 2c/ A S = 2(f^ ) 
which shows that £2 - f-] = 
and so on. So that for the n^^ null 
= nc/A S = nf-| 

and or more generally 




m - n AS 

where m >> n for best accuracy so that the measured 
difference in frequency between any two successive nulls can 
be related to the difference in path length required to 
produce that null by the equation 



^ _ C( receiver) 

AS 

where j is the actual sound velocity at 

the receiver at the time of reception. 
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IV. COMPARISON of the ANALYTICAL and COMPUTER SOLUTIONS for 
the ISOVELOCITY SOUND PROPAGATION CASE 

This analysis demonstrates the accuracy of the computer 
solution by using an isovelocity sound propagation model. 
The constant sound velocity condition results in straight 
line transmission paths as the sound travels from the source 
to the receiver. By solving the geometric model analytically 
and comparing the results to the computer output one can see 
the accuracy of the Romberg integration technique coupled 
with the four integrals derived in section III . The 
isovelocity case is the simplest propagation model yet it 
exercises all phases of the computer program's main 
calculation algorithm. The only section of the overall 
program not tested by the isovelocity case are the sound 
velocity calculations performed as the user inputs a varying 
temperature profile. Table IV-1 compares the difference in 
path lengths between the direct and reflected sound paths as 
calculated analytically and by the computer. Of the 20 data 
points compared, the computer missed the exact ray path 
length for a horizontal range of 7 50 meters by only 19.16 
centimeters in the worst case and averaged 5.02 centimeters 
error over all 20 cases. Table IV-2 shows the resulting 
horizontal range error due to the error of the computer's 
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calculation presented in Table IV-1. The average range error 
for the 20 points is only 1.52 meters. 

The computer program steps through sea depths of 0 to 
200 meters in 10 meter increments where each depth 
represents a possible target depth. To calculate a complete 
set of Snell constants, initial sound ray angles of 3 thru 
177 degrees in increments of 3 degrees { measured from the 
vertical ) are used. Ultimately the integration section 
derives a path length difference between the surface 
reflected ray and the direct path ray for each of these 
selected starting ray directions. The 59 initial rays times 
the 20 depth increments yields a total of 1,180 data points. 
These point solutions are used to interpolate values of path 
length difference for points every 10 meters of horizontal 
range out to 1000 meters which results in a two dimensional 
matrix containing 2000 ( 20 depth increments times 100 range 
increments ) horizontal range vs depth elements. 

To compare the computer solution to actual geometric 
results, five ranges for each of four depths were chosen at 
random. In all comparisons the receiver is set at a depth of 
27 meters (approximately 90 feet). 
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TABLE IV- 1 



SOURCE 

DEPTHS 

METERS 



10 



50 



90 



1 30 



170 



ISOVELOCITY CASE 
COMPARISON BETWEEN 
GEOMETRIC AND COMPUTER SOLUTIONS 



PATH LENGTH DIFFERENCE 
ANALYTICAL 
COMPUTER 



HORIZONTAL RANGES , METERS 



20 


100 


380 


750 


15.8107 


5. 1908 


1.4170 


0.71 95 


15.8228 


5.1939 


1.4418 


.6693 


49.0755 


23.5992 


7.0274 


3.5897 


49.0196 


23.6220 


7.0389 


3. 5714 


52.5987 


35.7218 


12.41 71 


6.4289 


52.6316 


35.71 43 


12.4688 


6.2267 


53.3450 


42.5840 


17.4439 


9.2169 


53.3808 


42 . 61 36 


17.5644 


9.0253 


53.6208 


46.431 2 


22 . 01 32 


11 .9302 


53.5714 


46.4396 


22.0264 


11 .7739 
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TABLE IV- 2 



COMPARISON of PROGRAM RANGES vs. ACTUAL RANGES 



ACTUAL 

SOURCE 



DEPTHS 

METERS 




HORIZONTAL 


RANGES, 


METERS 




20 


100 


380 


750 


10 


19.954 


99.935 


373.42 


806.3 


50 


20.134 


99.875 


379.37 


753. 87 


90 


19.745 


100.038 


378.33 


774. 84 


130 


19.42 


99.815 


377.075 


766.41 


170 


21 .26 


99.93 


379.725 


760.475 



COMPUTER 

HORIZONTAL RANGES, METERS 
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V. COMPUTER PROGRAM FLOW CHART 



A. FLOW CHART and PROGRAM BACKGROUND 

The following block diagram flow chart demonstrates the 
specific steps required to obtain an accurate horizontal 
range vs. depth solution using the Lloyd mirror interference 
phenomenon. Specific values or data used in the included 
program will be given at the bottom of each appropriate 
block. 

The actual numerical integration section of the program 
is taken directly from the Romberg integration program 
provided by Miller [Ref. 7: pp. 215-218]. 

The sound velocity equation used in the program was 
developed by Coppens [Refs. 6 and 8] and selected after 
several equations of varying complexity were considered. 
Specifically, equations derived by Wilson and Leroy [Ref. 1; 
p. 105], plus Mackenzie [Ref. 9] were compared with Coppens. 
All results were within 1 to 2 centimeters per second of 
each other. Coppen's equation was chosen because of the 
ability to correct for not only temperature, pressure and 
salinity, but also, latitude variations. 
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B . FLOW CHART 



1 . Initialization and Main Menu 
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DISPLAY MENU OF 
PROGRAM SUBROUTINES 
ALLOWING USER TO 
CHOOSE DESIRED 
OPERATION AS 
FOLLOWS : 




1 . 


REINITIALIZE RANDOM 
ACCESS ANSWER FILE 
AFTER INITIAL 
RUN AND PRIOR TO 
ALL SUBSEQUENT RUNS 


K!) 


^ • 


INPUT WATER TEMPERATURE 
INFO AS A FUNCTION OF 
DEPTH 


-0 


3. 


INPUT SEA SURFACE TEMP- 
ERATURE AND AVERAGE 
GRADIENT 


<D 


4. 


INPUT SOUND VELOCITY 
PROFILE AS A FUNCTION 
OF DEPTH 


*0 


5. 


COMPUTE HORIZONTAL RANGE 
VS. DEPTH SOLUTION AND 
STORE IN ANSWER FILE 





6. 


DISPLAY ANSWER FILE ON 
THE CRT SCREEN OR LIST 
ON A PRINTER 





7. 


EXIT PROGRAM AND RETURN 
TO OPERATING SYSTEM 


*<3 
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2. Reinitialize Random Access Answer File 
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3 . Input Water Temperature . Profile 
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4 



Input Sea Surface Temperature And Average Gradient 




(£) 
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5 



Input Sound Velocity Profile 



© 



INPUT SOUND VELOCITY 
IN METERS PER SECOND FOR 
INCREMENTS OF DEPTH IN METERS 
(PROGRAM USED 5 METER 
INCREMENTS) STORE VELOCITIES 
IN A TEMPORARY ARRAY 
(PROGRAM INTERPOLATES 1 METER 
INCREMENTS FROM 5 METER INPUTS) 
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. Com pute Range vs. Depth Solution And Store In Random 
Access File 




INPUT HYDROPHONE/RECEIVER 
DEPTH IN EITHER 
FEET OR METERS 



j LOOP THROUGH EACH POSSIBLE 
I SOURCE DEPTH, (THIS PROGRAM 
lUSES 10 METER INCREMENTS OF 
DEPTH DOWN TO 200 METERS) 

COMPUTING PATH LENGTH 
DIFFERENCES FOR EACH 
INCREMENT OF HORIZONTAL 
RANGE, (THIS PROGRAM USES 
10 METER INCREMENTS OF 
HORIZONTAL RANGE OUT 
TO 1000 METERS) 

BY DOING THE FOLLOWING: 

1 . DETERMINE THE UPPER 

AND LOWER LIMITS FOR THE 
FOUR SEPARATE NUMERICAL 
INTEGRATIONS 

2. LOOP THROUGH EACH 
POSSIBLE INITIAL RAY 

j ANGLE FROM EACH SOURCE 

I DEPTH COMPUTING THE FOLLOWING: 

a. RETRIEVE THE 

REQUIRED SOUND VELOCITIES 
FOR EACH INTEGRATION 

b . COMPUTE EACH 

UNIQUE SNELL CONSTANT 
GIVEN AN INITIAL RAY 
ANGLE AND SOUND VELOCITY 
FOR THE SOURCE DEPTH 

c. NUMERICALLY SOLVE 

THE FOUR INTEGRALS DERIVED 
IN THE MATHEMATICAL 
SOLUTION DISCUSSION, 

(cent. p. 9)- — 
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— — (cont. from p. 8) 

USING A ROMBERG IMPROVE- 
MENT [REF. 7] ON THE 
TRAPAZOIDAL NUMERICAL 
INTEGRATION TECHNIQUE 

d. IF INTEGRATION ERROR 
OCCURS WARN USER AND 
RETURN TO MENU 

e. STORE EACH INTEGRATION 
ANSWER IN THE APPROPRIATE 
TEMPORARY ARRAY 

f. CONTINUE TO LOOP THROUGH 
ALL POSSIBLE RAY ANGLES 
FROM THE SOURCE. (THIS 
PROGRAM USES ANGLES FROM 
3 DEGREES FROM THE VERT- 
ICAL TO 177 DEGREES FROM 
THE VERTICAL, IN 3 DEGREE 
INCREMENTS) 

3. USE LINEAR INTERPOLATION TO 
CALCULATE THE f ( DIFFERENCE ) 
FOR EACH INCREMENT OF RANGE 
FOR A GIVEN DEPTH 

4. STORE ANSWERS IN RANDOM 
ACCESS FILE 

5. CONTINUE TO LOOP THROUGH 
DEPTH INCREMENTS UNTIL 
DONE 



ON 



ERROR 




sir 

RETURN 
TO MENU 
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. Display Answer File On CRT Or List On Printer 



DISPLAY OUTPUT MENU TO 
USER ALLOWING FOR SELECTION 
OF: 

1 . CRT SCREEN 

2 . PRINTER 

RETRIEVE HORIZONTAL RANGE VS. 
DEPTH ANSWER MATRIX FROM 
RANDOM ACCESS FILE AND 
OUTPUT TO USERS CHOICE 
FROM ABOVE 




] RETURN 
I TO MENU 
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8 . Exit Program And Return To Operating System 
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• 

H 

> 


DEFINITION of COMPUTER PROGRAM VARIABLES 


ADD 


Equation used to compute raypath length in 
meters for both the direct and reflected paths. 


AFLAC 


Branching flag used in solving for answer matrix 
elements 


ANS 


Answer obtained from numeric integration 


AOUT 


Answer matrix element output to CRT screen or 
printer 


BFLAG 


Branching flag used in solving for answer matrix 
elements 


BNDLOW 


Lower limit for evaluating numeric integration 
II 


BOUNDHI 


Upper limit for evaluating numeric integration 


BOUNDLO 


Direct path lower limit 


CLAT 


Latitude variation used in sound velocity 
equation 


DDEEP 


First half of depth effect used in sound 
velocity equation 


DEEP 


Second half of depth effect used in sound 
velocity equation 


DEG. 


Average whole degree latitude of operating area 
used to calculate CLAT 


DELTA 


Line integral uniform divisions used in Romberg 
numeric integration routine 


DHI 


Direct path horizontal range upper bound used in 
interpolating direct path length for desired 
range increment 


DIRFRAC 


Direct path interpolation fraction 


DIRRAY 


Direct path interpolated path length 
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1 ^^ 











DLONG 

DLOW 

DPH 

DSHORT 

FANSWER 
FINAL ( 

GRAD 

HYDEEP 

HYDPTH 

HYSPD 

I 

IFLAG 

II 
IN 

INCR 

INN 

J 

JMP 

K 

M 

METER 



Dirct path length upper bound used in 
interpolating path length for desired range 
increment 

Direct path hoirzontal range lower bound used in 
interpolating direct path length for desired 
range increment 

Depth increment used for outputing answer matrix 

Direct path length lower bound used in 
interpolating path length for desired range 
increment 

FL1 Random access file used to store the horizontal 
range vs. depth answer matrix 

) Answer array used to temporarily hold the 
interpolated values before storing them in 
random access file 

Average temperature gradilent input by user 
Hydrophone/receiver depth in feet 
Hydrophone /receiver depth in meters 
Sound velocity at hydrophone/receiver depth 
Loop increment 

Loop increment used to branch through one of 
four integrated equations 

Loop increment 

Loop increment 

Loop increment 

Loop increment 

Counter 

Array element pointer 

Counter 

Counter 

Array element pointer 



46 






iin 



1 







MSTEP 

MSV 

MTRHI 

MTRLO 

MTRS 

MTRSTEP 

MULT 

N 

NANGL 

NMTR 

NN 

NSTEP 

NT 

NTRA 
NUM 
NUMB 
NUPLMT 
NX ( ) 

PART 

PIECES 

RAD 

RAYDIR ( 
RAYREF 



Loop increment 
Array element pointer 
Array element pointer 
Array element pointer 
Array element pointer 
Array element pointer 
Random access file pointer 
Array element pointer 
Ray angle from source 
Horizontal range increment 
Counter 

Ray angle from source 
Array element pointer 
Array element pointer 
Main menu index selection 
Output menu index selection 
Incremented source depth 

Array of counters used during integration 
accuracy testing loop 

Partitioning factor used in numerical 
integration routine 

Incremented divisor used to partition depth 
integration limits 

Conversion of ray angle from degrees to radians 

) Direct ray path length calculated from equation 
#12 page 25 

Reflected path length calculated from equation 
#13 page 25 
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REFFRAC 


Reflected path interpolation fraction 


REFLO 


Reflected path lower limit. Negative value of 
bound lo 


REFRAY 


Reflected path interpolated path length 


RHI 


Reflected path horizontal range upper bound used 
in interpolating reflected path length for 
desired range 


RLONG 


Reflected path length upper bound used in 
interpolating path length for desired range 
increment 


RLOW 


Reflected path horizontal range lower bound used 
in interpolating reflected path length for 
desired range 


RMTR 


Depth increment for interpolation routine 


RNGDIR ( ) 


Direct path horizontal range array calculated 
from equation §7 page 21 


RNGREF ( ) 


Reflected path horizontal range array calculated 
from equation #9 page 22 


RSHORT 


Reflected path length lower bound used in 
interpolating path length 


SAL IN 


Salinity value for operating area input by user 


SNDHI 


Sound velocity at numeric integral upper bound 


SNDLO 


Sound velocity at numeric integral lower bound 


SNDSPS ( ) 


Sound velocity profile array input directly by 
user or calculated from temperature, salinity, 
and latitude inputs 


SNELL 


Snell constant equation result calculated from 
initial ray angle and depth of source 


SOUND 


Sound velocity at incremented source depth 


SST 


Sea surface temperature 


START 


Initial estimate of integral answer 


SUM 


Equation used to compute horizontal range in 
meters for both the direct and reflected paths 
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SVP 


Input sound velocity profile by user 


TEMP 


Degrees Celsius temperature divided by 1 0 


TEST ( ) 


Successive numeric integration answer array used 
to check accuracy of answer 


TGTDPTH 


Incremented source depth 


TOL 


Accuracy desired for numerical integration 
routine answer 


TOT 


Summation of partitioned numerical ilntegration 
values 


TTEMP 


Input temperature used in calculating sound 
velocity profile 


WALK 


Integration routine smoothing factor 


X 


Incremented depth value for each partition 
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APPENDIX A 



Urick [Ref. 1: p. 124] briefly discusses the intensity of sound 
reaching a receiver via a direct path and a surface reflected path. The 
following derives the ratio of the total intensity of the signal received 
at a point when multipath Lloyd's mirror interference is present, to the 
total intensity resulting only from the direct path. 




WHERE : <I>, the time averaged intensity and I, the instantaneous 

intensity are defined as follows: 



(1) <i> 




T = 1 period 




( 2 ) 



(P„(r) + Pp(r))2 

2pc 



PqCt) and P[^(r) are both functions of time 
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GIVEN: 



’’received ^ ’’direct * ’’reflected 



0 

= ^ cos tot 
*^D 



«R 



cos [to(t-T)] 



where: u = surface pressure reflection coefficient 

T = time delay caused by the difference in the two path 
lengths 



( 2 ) 




P 



2 

0 



2pcR 



2 

D 



where: P = RMS pressure amplitude of the direct path ray 

NOTE : It is assumed that the transmission loss is due only to spherical 

spreading and that absorption in the medium is negligible. 
SOLUTION: 



P 

•^DIR 



0 , 
n COS \ COL ) 

% 



-uP 

’’ref " 



0 



COS[oo(t-T)] 



’’received ^ ’’direct ^ ’’reflected 



’’received ^ ’’o 



cos(o)t) - cos (oo(t-T)) 



2 P^ 

I = P_ = _0 

pC 2pC 



+ 4 C0S^(4t-x)) 



r2 
’^D 



- pPp cos (tot) cos [(o(t-T)] 

'^r'^d 
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SINCE; <I> = 4 /I dt 



AND 0) = 2Trf = |2. 






Ztt/o) 

/ 



T = ^ 

0 ) 



I dt 



<I> = 



0 ) P 



2irpC 



2tt/o) 

? / “ 



— / cos^(tot) dt + -^ 



<j 2Tr/a) 

I 



cos'^CooCt-t)) dt 



2u 



2it/(d 

/ ■ 



COS(o)t) COS(a)(t-x)) dt 



Solving the three integrals separately 



2ir/o) 



I cos‘‘(o)t) dt = [j + ^ sin(2 cot)] 



2tt/ ‘‘J 






IT 




52 



2 



2ir/ 

/• 



0 ) 



cos (w(t-t)) dt = 



_ M 



I + ^ sin( 2 a)(t-x)) 



2ir/(o 



_ U 






2 o)T)-sin(- 2 u,x)) 







2 ir/a) 

/ 



/ COS(o)t) COS(oj(t-T)) dt 
0 



_ 2 y 






2ir/ 

/ 



CO 



cos^(o)t) cos(toT)+sin(o)t) cos(u)t) sin(a)x) dt 



2 ir/a) 2 tt/o) 

= / cos^(»t) dt + 






i 



'^R^D 



/ 



sin(tot) cos(o,t) dt 



_ 2y cos(u)x) ,t . 1 
- R|,Rj| <2*5? sin(2»t) 



2 Tr/u) 



0 






'R‘'D 



2 ir/o) 

0 



_ 2u cos ( cot) /7T\ 2y sin((jOT) /p\ 
RrRd RrRq 



2ttu 



COS(o)T ) 
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mU 



Assembling the answers to the three integrals yields: 



U)P 



<I> = 



<I> = 



2irpC 

12_ 



<I> = 



o)Rd o)Rp R u 



V + ^ - IV COS(cut) 

Rn Rn '^r'^D 



u^Rp 2yRp 



2 pcR; 



NOTE: 



1 + — 2 — cos((^x) 



<Iq> = 



2pcR 



<I> 
Final ly 



<Iq> 



1 + 



2p2 

u Rq 



2 

D 

2pRr 



COS(o)t) 




This result shows the dependence of the intensity at the receiver on 
the difference in path lengths Rqjj^^qj and and the time delay, 

r, caused by this difference. This intensity ratio is a minimum (corre- 
sponding to a null on the Lloyd's mirror interference pattern) when 
cos COT = 1 or u)T = integer multiples of 2ir. 
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COMPUTER OUTPUT 



A sample of the actual computer output is listed below. 
Data values are frequency measurements Af (Hertz) between 
successive nulls on the Lloyd's mirror pattern (see Figure 
II-l) . 

Given a particular source depth and a Af, the source range 
can be determined from the corresponding output column. The 
source depth can be determined from the appropriate output row 
if the horizontal range to the source from the receiver and 
the Af are known. 



Group of Ranges 


From 10 


to 70 Meters 




Depth of 


10 Meters 










80.6 


94.8 114.1 


136.0 


159.8 


184.4 


210.1 


Depth of 


20 Meters 










41.9 


50.2 60.1 


71. 0 


82.7 


94.8 


107.3 


Depth of 


30 Meters 










31.6 


37.3 43.8 


50.8 


58. 3 


66.1 


74.3 


Depth of 


40 Meters 










29.2 


32.5 36.8 


41.7 


46.9 


52.5 


58. 3 


Depth of 


50 Meters 










28.5 


30.6 33.5 


36.9 


40.8 


44.9 


49.4 



55 



COMPUTER PROGRAM 



10 DEFINT I,J,K,L,M,N 

20 DIM NX(16) ,TEST( 1 36) ,SNDSPD( 201 ),RNGDIR( 180) 

30 DIM RNGREF(180),RAYDIR(180),RAYREF(180),FINAL(1 01 ) 

40 TOL = ,05 

50 OPEN "R", f^1 , "FANSWER.FL1'',4 
60 FIELD #1 , 4 AS FI $ 

100 PRINT CHR$(26): PRINT "SELECT OPERATION": PRINT 

no PRINT 1 /'INITIALIZE ANSWER FILE, DO AT THE BEGINNING OF 
EACH NEW BT INPUT" 

120 PRINT 2, "INPUT BT INFO IN 5 METER INCREMENTS" 

130 PRINT 3, "INPUT SEA SURFACE TEMP AND GRADIENT FOR FIRST 
200 METERS" 

135 PRINT 4, "INPUT SOUND VELOCITY PROFILE IN 5 METER 
INCREMENTS" 

140 PRINT 5, "CALCULATE LLOYD'S MIRROR HORIZONTAL RANGE ,VS. 
DEPTH CHART" 

150 PRINT 6, "OUTPUT CHART TO PRINTER" 

160 PRINT 7, "EXIT PROGRAM" 

170 PRINT:PRINT:INPUT "ENTER THE NUMBER OF THE OPERATION 
DESIRED ",NUM 

180 IF (NUM<1) OR (NUM>7) THEN PRINT "IMPROPER INPUT":G0T0 
100 

190 ON NUM GOSUB 1000,2000,5000,5500,6000,9600,10000 
200 GOTO 1 00 

1000 REM INITIALIZE ANSWER FILE 

1010 PRINT CHR${ 26):PRINT "ARE YOU SURE YOU WANT TO ZERO OUT 
THE ANSWER FILE" 

10 30 INPUT "ENTER Y FOR YES, N FOR NO ",A$ 

1040 IF A$o"Y" THEN RETURN 

1050 CLOSE #1 

1060 KILL "FANSWER.FL1" 

1070 OPEN "r",#1 ,"FANSWER.FL1",4 
1080 FIELD #1 , 4 AS FI $ 

1090 RETURN 

2000 REM INPUT BT INFO IN 5 METER INCREMENTS, THEN CALCULATE 
SOUND VELOCITY 

2010 REM USING ALAN R. COPPENS EQUATION PRESENTED IN J. 
ACOUS. SOC. OF AMER. 

2015 REM VOLUME 69, PAGE 862, 1981 AND PUBLISHED IN 

FUNDAMENTALS OF ACOUSTICS 

2016 REM THIRD EDITION, BY KINSLER FREY COPPENS AND 
SANDERS 

2020 GOSUB 3010: GOSUB 4000:PRINT CHR$(26) 

2030 PRINT "ENTER BT INFO IN 5 METER INCREMENTS STARTING 
WITH THE SURFACE" 

2040 FOR I = 0 TO 200 STEP 5 

2050 PRINT:PRINT I "M ETER S" : I NPUT "TEMP IN CELSIUS = 
",TTEMP:PRINT 

2060 PRINT "YOU INPUT "TTEMP" DEGREES CELSIUS": 
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INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 

";B$ 

2070 IF B$o"Y" THEN PRINT "REENTER TEMP":GOTO 2050 
2075 TEMP = TTEMP / 10! 

2080 METER = I 

2100 DEEP = (16.23 +.253*TEMP)* CLAT * CSNG( I) /1 000! +(.213- 

.1*TEMP) * 

(CLAT * (CSNG(I)/1 000!)) =*2 

2105 DDEEP = (.016 + .0002 * (SALIN-35)) * (SALIN- 

35)*TEMP*CLAT*CSNG(I)/1 000! 

2110 SNDSPD(METER) =1449.05 + 45.7 * TEMP - 5.21 * (TEMP'*2) + 

.23 * (TEMP^3) + 

(1 .333-.126*TEMP+.009*(TEMP^2) )*(SALIN- 
35)+(DEEP+DDEEP) 

2111 IF I = 0 GOTO 21 20 

2112 SLOW = SNDSPD( I - 5 ) 

21 1 3 SHI = SNDSPD( I ) 

21 1 4 FOR J = 1 TO 4 

21 15 SNDSPD( I - 5 + J ) = SLOW + ( J / 5 ) * ( SHI - SLOW ) 
21 1 6 NEXT J 
2120 NEXT I 
21 30 RETURN 

3000 REM ENTRY OF OPAREA LATITUDE 

3010 PRINT CHR$( 26): PRINT "ENTER AVERAGE DEGREES LATITUDE OF 
OPAREA" 

3020 INPUT "IN WHOLE DEGREES ",DEG 

3030 PRINT;PRINT "YOU ENTERED "DEG" DEGREES, IS THIS CORRECT 
?" 

3040 INPUT "ENTER Y FOR YES N FOR NO ",CC$ 

3050 IF CC$o"Y" THEN PRINT "REENTER LATITUDE ":GOTO 3010 
3060 CLAT = 1 - .0026 * (COS( DEG * 3.1 4 1 59 / 90 !)) 

3070 RETURN 

4000 REM ENTRY OF SALINITY VALUE FOR SOUND SPEED 

COMPUTATIONS 

4010 PRINT CHR$(26) 

4020 PRINT "ENTER THE SALINITY VALUE AS PARTS PER THOUSAND 
FOR THE OPAREA" 

4030 INPUT "IF UNKNOWN, ENTER A VALUE OF 35, SALINITY = 
" , SAL IN 

4040 PRINTzPRINT "THE SALINITY LEVEL IS "SALIN" PARTS PER 
THOUSAND" 

4050 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
",BB$ 

4060 IF BB$ <> "Y" GOTO 4010 
4070 RETURN 

5000 REM ENTER SEA SURFACE TEMP AND AVERAGE GRADIENT FOR 
FIRST 200 METERS 
5010 PRINT CHR$(26) 

5020 INPUT "ENTER THE SEA SURFACE TEMPERATURE IN DEGREES 
CELSIUS ",SST 

5030 PRINT:PRINT "YOU ENTERED "SST" DEGREES CELSIUS FOR THE 
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SEA SURFACE TEMP” 

50 40 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
”,C$ 

5050 IF C$ <> "Y" THEN PRINT "REENTER SEA SURFACE TEMP":GOTO 
5020 

5060 PRINT 

50 70 PRINT "enter THE TEMPERATURE GRADIENT (DEG 
CELSIUS /METER)" 

5075 INPUT "FOR THE FIRST 200 METERS ", GRAD: PRINT 

5080 PRINT- "YOU ENTERED A TEMPERATURE GRADIENT OF "GRAD" 

DEGREES CELSIUS /METER" 

5090 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
",D$ 

5100 IF D$ <> "Y" THEN PRINT "REENTER GRAD I ENT"; GOTO 5070 
5110 GOSUB 3010 
5120 GOSUB 4000 

5125 PRINT CHR$(26) ;PRINT "PLEASE WAIT, CALCULATING SOUND 
VELOCITY VALUES" 

5130 FOR I = 0 TO 200 

5140 TEMP = (SST + (CSNG(I) * GRAD)) / 101 

5150 DEEP = (1 6.23 + .253*TEMP)* CLAT * CSNG(I) /10001 + 

( .21 3-.1 *TEMP) 

* (CLAT * (CSNG(I) /1000! ))=*2 

5155 DDEEP = (.01 6 + .0002 * (SALIN - 35)) * (SALIN-35) * 
TEMP * CLAT * 

CSNG(I) / 1000! 

5170 SNDSPD(I) = 1449.05 + 45.7 * TEMP - 5.21 * (TEMP^2) + 
.23 * (TEMP=*3 ) 

+ ( 1 .333 - .1 26 * TEMP + . 0 0 9 *TE MP * 2 ) * ( S AL IN - 

35)+DEEP+DDEEP 
5180 NEXT I 

51 90 RETURN 

5510 REM ENTRY OF SOUND VELOCITY PROFILE 
5525 PRINT CHR$(26) 

5530 PRINT "ENTER SOUND VELOCITY PROFILE IN 5 METER 
INCREMENTS " 

5540 PRINT "STARTING WITH THE SURFACE SOUND VELOCITY IN 
METERS /SECOND" 

5550 FOR II = 0 TO 200 STEP 5 

5560 PRINT;PRINT II "METERS" ;INPUT " SOUND 
VELOCITY( METERS/ SEC) = ",SVP 

5570 PRINT; PRINT "YOU INPUT "SVP" METERS/SEC" 

5580 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
";SS$ 

55 90 IF SS$o"Y" THEN PRINT "REENTER SOUND VELOC ITY"; GOTO 
5560 

5610 SNDSPD(II) = SVP 

5611 IF II = 0 GOTO 5620 

5612 SLOW = SNDSPD( II - 5 ) 

5613 SHI = SNDSPD( II ) 

561 4 FOR JJ = 1 TO 4 
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5615 SNDSPD( II - 5 + JJ ) = SLOW + ( JJ / 5 ) * ( SHI - 

SLOW ) 

5616 NEXT JJ 
5620 NEXT II 
5630 RETURN 

6000 REM CALCULATE LLOYD'S MIRROR HORIZONTAL RANGE .VS. 
DEPTH CHART 

6005 PRINT:PRINT:PRINT "WILL HYDROPHONE DEPTH BE INPUT IN 
FEET OR METERS ?" 

6006 INPUT "ENTER F FOR FEET OR M FOR METERS ", AHYD$ 

6007 IF AHYD$ = "F" GOTO 6010 ELSE IF AHYD$ = "M" GOTO 6052 

6008 PRINT: PRINT: PRINT "INPUT ERROR" :GOTO 6005 

6010 PRINT CHR$(26) :INPUT "ENTER HYDROPHONE DEPTH IN FEET 
" ,HYDEEP 

6020 PRINT :PRINT "YOU ENTERED A HYDROPHONE DEPTH OF "HYDEEP" 
FEET" 

6030 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
",AA$ 

6040 IF AA$ <> "Y" GOTO 6010 

6050 HYDPTH = (HYDEEP * .3048) 

6051 GOTO 6060 

6052 PRINT CHR$( 26): INPUT "ENTER HYDROPHONE DEPTH IN METERS 
", HYDPTH 

6053 PRINT:PRINT "YOU ENTERED A HYDROPHONE DEPTH OF 
"HYDPTH" METERS" 

60 54 INPUT " IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
",BHY$ 

6055 IF BHY$ <> "Y" GOTO 6052 
6060 HYSPD = SNDSPD( HYDPTH ) 

61 70 FOR NUPLMT = 10 TO 200 STEP 1 0 
6175 PRINT "NUPLMT=" NUPLMT 

6180 TGTDPTH = CSNG( NUPLMT ) 

6185 IF HYDPTH = TGTDPTH THEN HYDPTH = HYDPTH - 3 

6190 IF TGTDPTH > HYDPTH GOTO 6230 

6200 BOUNDLO = TGTDPTH 

6210 BOUNDHI = HYDPTH 

6220 GOTO 6250 

6230 BOUNDLO = HYDPTH 

6240 BOUNDHI = TGTDPTH 

6250 REFLO = -(BOUNDLO) 

6360 FOR NSTEP = 3 TO 1 77 STEP 3 

6365 PRINT "NSTEP=" NSTEP 

6370 RAD = CSNG(NSTEP) * 3.1 41 59 / 1 801 

6380 MTRSTEP = TGTDPTH 

6390 SOUND = SNDSPD( MTRSTEP) 

6400 MTRHI = BOUNDHI 
6410 SNDHI = SNDSPD( MTRHI) 

6420 MTRLO = BOUNDLO 
6430 SNDLO = SNDSPD(MTRLO) 

64 40 SNELL = SIN( RAD ) / SOUND 
6480 FOR IFLAG = 1 TO 4 
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7000 REM ROMBERG NUMERICAL INTEGRATION SUBROUTINE 
7015 START = 01 
7020 PIECES = II 
7030 NX(1 ) = 1 

7040 IF (IFLAG=2) OR (IFLAG=4) THEN BNDLOW=REFLO ELSE 
BNDLOW=BOUNDLO 

7041 IF (IFLAG = 2) OR (IFLAG = 4) GOTO 7045 

7042 IF (TGTDPTH > HYDPTH) AND (NSTEP > 90) GOTO 8105 

7043 IF (TGTDPTH < HYDPTH) AND (NSTEP < 90) GOTO 8105 

7044 GOTO 7050 

7045 IF NSTEP > 90 GOTO 8105 

7050 DELTA = (BOUNDHI - BNDLOW) / PIECES 

7051 DENOMl = (SNELL) *2 * (SNDLO)*2 

7052 DENOM2 = (SNELL)*2 * (SNDHI)^2 

7053 IF DENOMl = II THEN DENOMl = .999 ;IF DENOM2 = II THEN 
DENOM2 = .999 

7054 IF ( DENOMl > II ) OR ( DENOM2 > II ) GOTO 8105 
7060 IF (IFLAG=1) OR (IFLAG=2) GOTO 7085 

7065 FOR IN = 1 TO 2 

70 68 IF IN = 1 THEN DENOM = DENOMl ELSE DENOM = DENOM2 
7070 START = START + ((II)/ SQR(1 I - DENOM)) 

7075 NEXT IN 

7080 GOTO 7100 

7085 FOR INN = 1 TO 2 

7088 IF INN = 1 THEN DENOM = DENOMl ELSE DENOM = DENOM2 

7089 IF INN = 1 THEN SPEED = SNDLO ELSE SPEED = SNDHI 

7090 START = START + ((SNELL * SPEED) / SQR(1 I - DENOM)) 

7095 NEXT INN 

71 00 TOT = START / 2 1 

7110 TEST(l) = DELTA * START / 21 

7120 N = 1 

71 30 NN = 2 

7140 N = N + 1 

7150 WALK = 41 

71 60 NX(N) = NN 

7170 MSTEP = CINT( PIECES ) 

7180 PIECES = PIECES * 21 

7190 DELTA = (BOUNDHI - BNDLOW) / PIECES 

7500 REM COMPUTE TRAPEZOIDAL SUM 

7510 FOR INCR = 1 TO MSTEP 

7520 PART = CSNG(INCR) * 21 - II 

7530 X = BNDLOW + (DELTA * PART) 

7540 MTRS = CINT ( ABS (X) ) 

7545 IF MTRS > 200 THEN MTRS = 200 
7550 SPEED = SNDSPD( MTRS ) 

7552 DENOM3 = (SNELL)*2 * (SPEED) *2 

7554 IF DENOM3 = II THEN DENOM3 = .999 

7555 IF DENOM3 > II GOTO 8105 

7570 IF (IFLAG=1) OR (IFLAG=2) GOTO 7600 
7580 TOT = TOT + ( (II) / SQR(1I - DENOM3)) 

7590 GOTO 7610 



60 





iflffltlWi 



DEN0M3)) 



7600 TOT = TOT + ( (SNELL * SPEED) / SQR{1! - 
7610 NEXT INCR 

76 20 TEST{ NN ) = TOT * DELTA 
7630 K = N - 1 
76 40 NTRA = NX{ K ) 

7650 FOR M = 1 TO K 
76 60 J = NN + M 
76 70 NT = NTRA + M - 1 
7680 TEST(J) = ( TEST(J-I) * WALK - TEST(NT)) / (WALK-1!) 

76 90 WALK = WALK * 41 
7700 NEXT M 

7705 IF (IFLAG = 1 OR IFLAG = 2) AND (TEST(J) > 3000!) GOTO 
81 05 

7710 IF N < 5 GOTO 7760 

7720 IF TEST( NN + 1 ) = 01 GOTO 7750 

7730 IF ABS(TEST(NTRA+ 1 )-TEST(NN + 1)) <= (ABS(TEST(NN + 

1)*TOL)) GOTO 8000 

7740 IF ABS( TEST(NN-I) - TEST(J)) <= (ABS( TEST(J) * TOL)) 
GOTO 8000 

7750 IF N > 15 GOTO 9500 
7760 NN = J + 1 
7770 GOTO 7140 
8000 ANS = TEST(J) 

8030 ON IFLAG GOTO 8040,8060,8080,8100 

8040 RNGDIR(NSTEP) = ANS 

8050 GOTO 8105 

8060 RNGREF(NSTEP) = ANS 

8070 GOTO 8105 

8080 RAYDIR(NSTEP) = ANS 

8090 GOTO 8105 

8100 RAYREF(NSTEP) = ANS 

81 05 NEXT IFLAG 

8110 NEXT NSTEP 

8115 REM USE LINEAR INTERPOLATION TO CALCULATE ANSWER CHART 
VALUES IN RANGE 

8117 REM AND DEPTH INCREMENTS OF 10 METERS 

8120 FOR NMTR = 1 TO 100 

8130 RMTR = CSNG( NMTR ) * 10! 

8135 DIRRAY = 01 ; REFRAY = 0! 

81 40 AFLAG = 0! 

81 50 BFLAG = Ol 

8160 IF TGTDPTH < HYDPTH GOTO 8179 

8161 FOR NANGL = 3 TO 87 STEP 3 

8162 IF RNGDIR( NANGL ) <= RMTR GOTO 8171 

8163 DLOW = RNGDIR( NANGL - 3 ):DHI = RNGDIR( NANGL ) 

8164 IF DLOW = DHI THEN DIRFRAC = 0 ;GOTO 8166 

8165 DIRFRAC = (RMTR - DLOW) / (DHI - DLOW) 

81 66 IF DLOW > 0 GOTO 81 68 

8167 DSHORT = BOUNDHI - BOUNDLO :GOTO 8169 

8168 DSHORT = RAYDIR( NANGL - 3 ) 

81 69 DLONG = RAYDIR( NANGL ) 
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81 70 DIRRAY = DSHORT + ( DLONG - DSHORT ) * DIRFRAC : GOTO 
821 0 

8171 IF ( NANGL = 87 ) AND ( RNGDIR( 87 ) > 0 ) GOTO 8172 
ELSE GOTO 8177 

8172 DLOW = RNGDIR( 87 ) 

8173 IF DLOW = RMTR THEN DIRRAY = RAYDIR( 87 ) : GOTO 8210 

8174 DIRFRAC = ( RMTR - DLOW ) / ( 1000 - DLOW ) 

8175 DSHORT = RAYDIR( NANGL ) : CLONG = SQR ( 10 00 "*2 + 

(BOUNDHI - BOUNDLO)=2 ) 

8176 DIRRAY = DSHORT + ( CLONG - DSHORT ) * DIRFRAC : GOTO 
8210 

8177 NEXT NANGL 

8178 GOTO 8410 

8179 FOR MANGL = 93 TO 177 STEP 3 

8180 IF RMTR >= RNGDIR( MANGL ) GOTO 8181 ELSE GOTO 8195 

8181 IF ( MANGL = 93 ) AND ( RNGDIR( 93 ) > 0 ) GOTO 8182 
ELSE GOTO 8186 

8182 DLOW = RNGDIR( 93 ) 

8183 IF DLOW = RMTR THEN DIRRAY = RAYDIR( 93 ) ; GOTO 8210 

8184 DIRFRAC = ( RMTR - DLOW ) / ( 1 000 - DLOW ) ; DSHORT = 

RAYDIR( 93 ) : 

BLONG = SQR ( 1 00 0=2 + ( BOUNDHI - BOUNDLO ) = 2 ) 

8185 DIRRAY= DSHORT + ( BLONG - DSHORT ) * DIRFRAC : GOTO 

8210 

8186 DLOW = RNGDIR( MANGL ) :DHI = RNGDIR( MANGL - 3 ) 

8187 IF ( DLOW <= RMTR ) AND ( DHI >= RMTR ) GOTO 8188 ELSE 
GOTO 81 95 

8188 IF DLOW = DHI THEN DIRFRAC = 0 : GOTO 8190 

8189 DIRFRAC = ( RMTR - DLOW ) / ( DHI - DLOW ) 

81 90 IF DLOW > 0 GOTO 81 92 

8191 DSHORT = BOUNDHI - BOUNDLO : GOTO 8193 

8192 DSHORT = RAYDIR( MANGL ): DLONG = RAYDIR( MANGL - 3) 

81 94 DIRRAY = DSHORT + ( DLONG - DSHORT ) * DIRFRAC : GOTO 
8210 

8195 IF ( MANGL <> 177 ) GOTO 8199 

8196 DLOW = 0 ; DHI = RNGDIR( 1 77 ) : DIRFRAC = RMTR / DHI 

8197 DSHORT = BOUNDHI - BOUNDLO : DLONG = RAYDIR( 177 ) 

81 98 DIRRAY = DSHORT + ( DLONG - DSHORT ) * DIRFRAC : GOTO 

821 0 

8199 NEXT MANGL 

8200 GOTO 8410 

8205 IF DHI = DLOW THEN DIRFRAC = 0 : GOTO 8211 

8210 FOR NREF = 3 TO 87 STEP 3 

8211 IF RNGREF( NREF ) <= RMTR GOTO 8220 

8212 RLOW = RNGREF( NREF - 3 ) ; RHI = RNGREF( NREF ) 

8213 IF RLOW = RHI THEN REFFRAC = 0 : GOTO 8215 

8214 REFFRAC = ( RMTR - RLOW ) / ( RHI - RLOW ) 

8215 IF RLOW > 0 GOTO 8217 

8216 RSHORT = BOUNDHI + BOUNDLO : GOTO 8218 

8217 RSHORT = RAYREF( NREF - 3 ): RLONG = RAYREF(NREF) 

8219 REFRAY = RSHORT + ( RLONG - RSHORT ) * REFFRAC :GOTO 
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8400 

8220 IF ( NREF = 87 ) AND (RNGREF(87) > 0 ) GOTO 82 21 ELSE 
GOTO 8226 

8221 RLOW = RNGREF( 87 ) 

8222 IF RLOW = RMTR THEN REFRAY = RAYREF{ 87 ): GOTO 8400 

8223 REFFRAC = ( RMTR - RLOW ) / ( 1 000 - RLOW ) 

8224 RSHORT = RAYREF( NREF ) : ALONG = SQR ( 1 000'*2 + ( 
BOUNDHI + BOUNDLO )^2 ) 

82 25 REFRAY = RSHORT + ( ALONG - RSHORT ) * REFFRAC : GOTO 

8400 

8226 NEXT NREF 

8227 GOTO 8410 

83 70 REM OUTPUT FINAL ANSWERS FOR EACH DEPTH AND RANGE TO 
DISK FILE 

8385 LPRINT "HYSPD=" HYSPD 

8400 FINAL( NMTR ) = HYSPD / ( REFRAY - DIRRAY ) 

8410 NEXT NMTR 

8420 MULT = (TGTDPTH ±10) - 1 

8430 FOR JMP = 1 TO 100 

8440 LSET FI $ = MKS$( FINAL( JMP) ) 

8450 RECNUM% = JMP + (MULT * 1 00 ) 

8460 PUT RECNUM% 

8470 NEXT JMP 

8480 FOR NUL = 0 TO 100 

8490 FINAL(NUL) = 01 

8500 NEXT NUL 

8510 FOR NUL = 0 TO 135 

8520 TEST(NUL) = 01 

8530 NEXT NUL 

8540 FOR NUL = 0 TO 1 5 

8550 NX(1 5) = 01 

8560 NEXT NUL 

8570 FOR NUL = 0 TO 1 79 

85 80RNGDIR(NUL)=0!:RNGREF(NUL)=0:RAYDIR(NUL) = 0: 

RAYREF ( NUL ) =0 
8620 NEXT NUL 
9000 NEXT NUPLMT 
9010 RETURN 

9500 PRINT; PRINT: PRINT "INPUT DATA HAS CAUSED AN ERROR" 

9510 GOTO 100 

9600 REM PRINT OUT GRAPH OF HORIZONTAL RANGE .vs. DEPTH 
9610 PRINT CHR$(26) 

9620 PRINT ;PRINT;PRINT " SELECT OUTPUT DEVICE ": PRINT 

9630 PRINT 1," VIDIO SCREEN PRINT 
9640 PRINT 2," PRINTER PRINT 

9650 INPUT " ENTER THE NUMBER FOR OUTPUT DEVICE DESIRED 
" ,NUMB 

96 60 IF ( NUMB < 1 OR NUMB > 2 ) THEN PRINT "IMPROPER INPUT 
";GOTO 9620 
9665 PRINT; PRINT 

96 70 PRINT " THE OUTPUT MATRIX IS PRESENTED IN 10 METER 
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INCREMENTS FOR BOTH " 

9680 PRINT " HORIZONTAL RANGES AND FOR DEPTHS OF 10 METERS 
TO 200 METERS. " 

96 90 PRINT " HORIZONTAL RANGES FROM THE RECEIVER TO THE 
TARGET , IN " 

9700 PRINT " GROUPS OF 100 METERS, INCREASE FROM LEFT TO 
RIGHT, " 

9710 PRINT '' ie. 10 TO 100 METER RANGES WILL BE DISPLAYED 
FIRST, FOLLOWED" 

9720 PRINT "BY 110 TO 200 METERS AND SO ON UP TO 1000 
METERS . " 

9730 PRINT " DEPTHS OF 0 TO 200 METERS WILL BE DISPLAYED 
FOR EACH GROUP OF " 

9740 PRINT " RANGES, WITH THE FIRST ROW REPRESENTING 10 
METERS, THE " 

9750 PRINT " NEXT ROW 20 METERS, ON DOWN THE PAGE UNTIL ROW 
20 IS " 

9760 PRINT " PRINTED FOR A DEPTH OF 200 METERS."; PRINT: PRINT 

9770 FOR I = 0 TO 9 

9780 THIS = 10 + 100 * I 

9790 THAT = 1 00 + 1 00 * I 

9800 IF NUMB = 1 THEN PRINT; 

PRINT "GROUP OF RANGES FROM "THIS" TO 
"THAT"METERS";GOTO 9 82 0 
9805 LPRINT;LPRINT; 

9810 LPRINT;LPRINT "GROUP OF RANGES FROM "THIS" TO "THAT" 
METERS" 

9820 FOR J = 0 TO 19 
9830 DPH = 10 + 10 * J 

9840 IF NUMB = 1 THEN PRINT:PRINT "DEPTH OF "DPH" 

METERS" .‘GOTO 9860 

9850 LPRINT;LPRINT "DEPTH OF "DPH" METERS" 

9860 FOR K = 1 TO 10 

9870 RECNUM% = K + 100 * J + 10 * I 

9880 GET #1,RECNUM% 

98 90 AOUT = CVS{ FI $ ) 

9900 IF NUMB = 2 GOTO 9930 

9910 PRINT USING ";AOUT; 

9920 GOTO 9940 

9930 LPRINT USING ";AOUT; 

9940 NEXT K 
9950 NEXT J 
9960 NEXT I 
9965 LPRINT 
9970 RETURN 

10000 REM CLOSE ALL FILES AND EXIT PROGRAM 

10010 PRINT CHR$( 26); PRINT "DO YOU REALLY WANT TO QUIT" 

10 020 INPUT "ENTER A Y FOR YES, AN N FOR NO ",Z$ 

10030 IF Z$ = "Y" THEN END ELSE RETURN 
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